
clear

bootidx = [1:80];

SaveCon = 'Results/full_CF.mat';
load(SaveCon);
CFest = CF;
SaveRes = 'Results/full.mat';
load(SaveRes);

% Recover FID_Grid
optCount = opt; optCount.SelSubSample = 0; 
optCount.TrunkAge = 0; % Option to trucate age (added later, so may not be present in some results files)
[FID_Grid,Action,Age,AgeState,States,pproc,DS] = prepdataM10(optCount);

% Make sure FID_Grid match:
GridNMatch = any(CF.FID_Grid - FID_Grid ~= 0);
if GridNMatch
    disp(' CF grid does not match with estimates option generated grid ')
end
% Load Bootstrap elasticities estimates:
AcreageSweight_b = nan(size(bootidx)); AcreageAweight_b = nan(size(bootidx)); YieldSweight_b = nan(size(bootidx)); 
Output_b = nan(size(bootidx)); YieldShare_b = nan(size(bootidx));
for b = bootidx
    b
    cfbootname = ['Bootstrap/CF/CF_full_b', num2str(b), '.mat'];
    if exist(cfbootname) ~= 2; continue; end
    load(cfbootname)
    IndAcreageElas = ((CF.LandUseCF05./CF.LandUseBase) - 1)/0.05;
    IndYieldElas   = (((CF.OutputDetailCF05./CF.LandUseCF05)./(CF.OutputDetailBase./CF.LandUseBase))-1)/0.05;
    IndOutputElas  = ((CF.OutputDetailCF05./CF.OutputDetailBase) - 1)/0.05;
    
    
    AcreageSweight_b(b)   = [nansum(IndAcreageElas.*CF.OutputDetailBase)/nansum(CF.OutputDetailBase)];
    AcreageAweight_b(b)   = [nansum(IndAcreageElas.*CF.LandUseBase)/nansum(CF.LandUseBase)];
    YieldSweight_b(b)     = [nansum(IndYieldElas.*CF.OutputDetailBase)/nansum(CF.OutputDetailBase)];
    Output_b(b)           = [nansum(IndOutputElas.*CF.OutputDetailBase)/nansum(CF.OutputDetailBase)];
    YieldShare_b(b)       = [YieldSweight_b(b)/Output_b(b)];
end
%% Elasticities decomposition

% Individual acreage and yield elasticities:
IndAcreageElas = ((CFest.LandUseCF05./CFest.LandUseBase) - 1)/0.05;
IndYieldElas   = (((CFest.OutputDetailCF05./CFest.LandUseCF05)./(CFest.OutputDetailBase./CFest.LandUseBase))-1)/0.05;
IndOutputElas  = ((CFest.OutputDetailCF05./CFest.OutputDetailBase) - 1)/0.05;

TableElas.AcreageSweight = [nansum(IndAcreageElas.*CFest.OutputDetailBase)/nansum(CFest.OutputDetailBase); nanstd(AcreageSweight_b)];
TableElas.AcreageAweight = [nansum(IndAcreageElas.*CFest.LandUseBase)/nansum(CFest.LandUseBase); nanstd(AcreageAweight_b)];
TableElas.YieldSweight   = [nansum(IndYieldElas.*CFest.OutputDetailBase)/nansum(CFest.OutputDetailBase);nanstd(YieldSweight_b)];
TableElas.Output         = [nansum(IndOutputElas.*CFest.OutputDetailBase)/nansum(CFest.OutputDetailBase);nanstd(Output_b)];
TableElas.YieldShare     = [TableElas.YieldSweight(1)/TableElas.Output(1); nanstd(YieldShare_b)];

TableElas = struct2table(TableElas);
TableElas.Properties.RowNames = {'Estimate', 'B. Std. Err.'}


%% Graph short-medium run elasticities

Time = 1:150;
for t = Time
    Output(t)  = (mean(CFest.OutputPathCF05(t,:),2)/mean(CFest.OutputPathBase(t,:),2) - 1)/0.05;
    Yield(t)   = (mean(CFest.YieldPathCF05(t,:),2)/mean(CFest.YieldPathBase(t,:),2) - 1)/0.05;
    Acreage(t) = (mean(CFest.AcreagePathCF05(t,:),2)/mean(CFest.AcreagePathBase(t,:),2) - 1)/0.05;
end

figure1 = figure(1); % Short and long run
% Create axes
axes1 = axes('Parent',figure1,'XTick',[1 3 5 7 9 11 13 15 20 25 30 40 50 60 70 80 90]);
xlim(axes1,[1 90]);
box(axes1,'on');
grid(axes1,'on');
hold(axes1,'all');

% Create multiple lines using matrix input to plot
plot1 = plot([Output(1:90)' Yield(1:90)' Acreage(1:90)'],'Parent',axes1,'LineWidth',3);
set(plot1(1),'DisplayName','Output');
set(plot1(2),'DisplayName','Yield');
set(plot1(3),'DisplayName','Acreage');

% Create xlabel
xlabel('years');

% Create ylabel
ylabel('elasticity');

% Create legend
legend1 = legend(axes1,'show');

figure2 = figure(2); % Short run
% Create axes
axes1 = axes('Parent',figure2,'XTick',[1 3 5 7 9 11 13 15 17]);
xlim(axes1,[1 17]);
box(axes1,'on');
grid(axes1,'on');
hold(axes1,'all');

% Create multiple lines using matrix input to plot
plot1 = plot([Output(1:17)' Yield(1:17)' Acreage(1:17)'],'Parent',axes1,'LineWidth',3);
set(plot1(1),'DisplayName','Output','Color','black','LineStyle','-');
set(plot1(2),'DisplayName','Yield','Color','black','LineStyle','--');
set(plot1(3),'DisplayName','Acreage','Color','black','LineStyle',':');

% Create xlabel
xlabel('years');

% Create ylabel
ylabel('elasticity');

% Create legend
legend2 = legend(axes1,'show');


%% Acreage Elasticity Decomposition:

WildVariable = {'ECO_NUM'};
WildIdx = double(DS(:,WildVariable));
cf = 1; % Counterfactual for elasticities.
TableWildElasticity.WildIdx        = unique(WildIdx);
TableWildElasticity.Elasticity     = zeros(length(unique(WildIdx)),1);
TableWildElasticity.BaselineShare  = zeros(length(unique(WildIdx)),1);
TableWildElasticity.WildTotal      = zeros(length(unique(WildIdx)),1);
TableWildElasticity.TotElasContrib = zeros(length(unique(WildIdx)),1);
TableWildElasticity.ShareOfExpan   = zeros(length(unique(WildIdx)),1);

w_count = 1;
for w = unique(WildIdx)'
    TableWildElasticity.Elasticity(w_count)     = nansum(IndAcreageElas(WildIdx == w).*CFest.LandUseBase(WildIdx == w))/nansum(CFest.LandUseBase(WildIdx == w));
    TableWildElasticity.TotElasContrib(w_count) = nansum(IndAcreageElas(WildIdx == w).*CFest.LandUseBase(WildIdx == w))/nansum(CFest.LandUseBase);
    TableWildElasticity.ShareOfExpan(w_count)   = TableWildElasticity.TotElasContrib(w_count)/TableElas.AcreageAweight(1);
    TableWildElasticity.BaselineShare(w_count)  = nanmean(CFest.LandUseBase(WildIdx == w));
    TableWildElasticity.WildTotal(w_count)      = sum(WildIdx == w);
    w_count = w_count + 1;
end

TableWildElasticity = struct2table(TableWildElasticity)

%% Cropland and Pasture effects of price changes:

Pasture  = DS.nasapasture2000;
PastureLandUse = Pasture.*CFest.LandUseBase;
Crop = DS.nasacropland2000;
CropLandUse = Crop.*CFest.LandUseBase;

TableCroplandPasture.Cover = {'Pasture';'Cropland'};
TableCroplandPasture.Elasticity = zeros(2,1);
TableCroplandPasture.TotElasContrib = zeros(2,1);
TableCroplandPasture.ShareOfExpan = zeros(2,1);
TableCroplandPasture.BaselineShare = zeros(2,1);
TableCroplandPasture.TotalArea     = zeros(2,1);

TableCroplandPasture.Elasticity(1) = nansum(IndAcreageElas.*PastureLandUse)/nansum(PastureLandUse);
TableCroplandPasture.Elasticity(2) = nansum(IndAcreageElas.*CropLandUse)/nansum(CropLandUse);
TableCroplandPasture.TotElasContrib(1) = nansum(IndAcreageElas.*PastureLandUse)/nansum(CFest.LandUseBase);
TableCroplandPasture.TotElasContrib(2) = nansum(IndAcreageElas.*CropLandUse)/nansum(CFest.LandUseBase);
TableCroplandPasture.ShareOfExpan(1) = TableCroplandPasture.TotElasContrib(1)/TableElas.AcreageAweight(1);
TableCroplandPasture.ShareOfExpan(2) = TableCroplandPasture.TotElasContrib(2)/TableElas.AcreageAweight(1);
TableCroplandPasture.BaselineShare(1) = nansum(PastureLandUse)/nansum(Pasture);
TableCroplandPasture.BaselineShare(2) = nansum(CropLandUse)/nansum(Crop);
TableCroplandPasture.TotalArea(1)     = nansum(Pasture);
TableCroplandPasture.TotalArea(2)     = nansum(Crop);

TableCroplandPasture = struct2table(TableCroplandPasture)

%% Marginal fields (how new fields compare to old fields)

NewFields = CFest.LandUseCF05-CFest.LandUseBase;
NewFields(NewFields < 0) = 0;

PotYield = double(DS(:,opt.PotenYield));
PotYieldNew = sum(NewFields.*PotYield)/sum(NewFields);
PotYieldOld = sum(CFest.LandUseBase.*PotYield)/sum(CFest.LandUseBase);

TransCost = double(DS(:,opt.TransCostVar));
TransCostNew = sum(NewFields.*TransCost)/sum(NewFields);
TransCostOld = sum(CFest.LandUseBase.*TransCost)/sum(CFest.LandUseBase);



